drop _all
capture program drop SPR


program define SPR, rclass
{
  syntax [, n(real 50000)]
  tempvar e u x xd d y0 y1  ys c y del one lny lnyhat res2 estw250 estw300 estw350 d_i estlc025 estlc030 estlc035 estec025 estec030 estec035
  
set obs 40000

gen `e'=rnormal() 
gen `x'=2*runiform() 
gen `d'=(`x'+`e')>1
gen `xd'=`x'*`d'
gen `u'=rexponential(1)
scalar betd=0.5
gen `y0'=`x'+`u'
gen `y1'=`y0'+betd
gen `ys'=`y0'+(`y1'-`y0')*`d'
gen `c'=rexponential(8)
gen `y'=`ys'*(`ys'<`c')+`c'*(`c'<`ys')
gen `one'=1
gen `del'=(`ys'<`c')


// true effect computed at t evaluation points (t=2.5 3 3.5)
forv i=25(5)35 {
**# Bookmark #1
	scalar id=`i'/10
	qui gen true0=exp(-(id-`x'))
	qui gen true1=exp(-(id-betd-`x'))
    qui gen mtrue=(true1/true0)-1
	qui sum mtrue
    scalar true_`i'= r(mean)
   drop true0 true1 mtrue
}


forv i=25(5)35 {

	scalar id=`i'/10
	
	gen d_i=(`y'>id)
    qui reg d_i `x' `d' `xd'
    scalar b00=_b[_cons]
	scalar b10=_b[`x']
	scalar b20=_b[`d']
	scalar b30=_b[`xd']
	scalar eps=0.00001

	
//********** Estimator specifying linear cumulative hazard at each t **************//	
	mlexp ( (1-d_i)*ln(1-(exp(-(({b0=b00}+{b1=b10}*`x'+{b2=b20}*`d'+{b3=b30}*`xd')>eps)*({b0=b00}+{b1=b10}*`x'+{b2=b20}*`d'+{b3=b30}*`xd')+(({b0=b00}+{b1=b10}*`x'+{b2=b20}*`d'+{b3=b30}*`xd')<eps)*eps))) -d_i*abs({b0=b00}+{b1=b10}*`x'+{b2=b20}*`d'+{b3=b30}*`xd')  )

	matrix blc=e(b)
	scalar n=e(N) 
	scalar blc0=_b[/b0]
	scalar blc1=_b[/b1]
	scalar blc2=_b[/b2]
	scalar blc3=_b[/b3]
    gen `estlc0`i''=exp(-(blc2+blc3*`x'))-1 
	qui sum `estlc0`i''
	scalar estlc`i'=r(mean)
	
//************** Estimator specifying exponential cumulative hazard at each t **************//	
	mlexp ( d_i*ln(exp(-exp({b0=b00}+{b1=b10}*`x'+{b2=b20}*`d'+{b3=b30}*`xd'))) + (1-d_i)*ln(1-(exp(-exp({b0=b00}+{b1=b10}*`x'+{b2=b20}*`d'+{b3=b30}*`xd'))))  )

	matrix bec=e(b)
	scalar bec0=_b[/b0]
	scalar bec1=_b[/b1]
	scalar bec2=_b[/b2]
	scalar bec3=_b[/b3]

	gen `estec0`i''= exp(-exp(bec0+bec1*`x'+bec2+bec3*`x')+exp(bec0+bec1*`x')) -1 
	qui sum `estec0`i''
	scalar estec`i'=r(mean)
	
	
	mata {
	   st_view(y=.,.,tokens("d_i"))
	   st_view(mx=.,.,tokens("`one' `x' `d' `xd'"))
	   st_view(m=.,.,tokens("`one' `x'"))
	   st_view(md=.,.,tokens("`d' `xd'"))  
	   st_view(mint=.,.,tokens("`one' `x'")) // x list used as interaction term
	   st_view(one=.,.,tokens("`one' "))
	   
	   blc=st_matrix("blc")
	   n=rows(y)
	   bmdlc=(blc[1,3]\blc[1,4])
	   slc=-(y:*mx):/(one-exp(-mx*blc'))
	   etalc=slc*invsym(slc'*slc/n)
	   gglc=exp(-mint*bmdlc)
	   el1lc=(0\0) 
	   el2lc=mean(-mint:*gglc)
	   ellc=(el1lc\el2lc')
	   ggmlc=mean(gglc)
	   inflc=gglc-one*(ggmlc)+etalc*ellc
	   vlc=inflc'*inflc/(n^2)
	   st_matrix("vlc", vlc)
	   
	   bec=st_matrix("bec")
	   bmdec=(bec[1,3]\bec[1,4])
	   bmec=(bec[1,1]\bec[1,2])
	   sec=-(y:*exp(mx*bec'):*mx):/(one-exp(-exp(mx*bec')))
	   etaec=sec*invsym(sec'*sec/n)
	   ggec=exp( -exp(m*bmec+mint*bmdec)+exp(m*bmec))
	   el1ec=mean(-m:*ggec:*(exp(m*bmec+mint*bmdec)-exp(m*bmec))) 
	   el2ec=mean(-mint:*ggec:*exp(m*bmec+mint*bmdec))
	   elec=(el1ec'\el2ec')
	   ggmec=mean(ggec)
	   infec=ggec-one*(ggmec)+etaec*elec
	   vec=infec'*infec/(n^2)
	   st_matrix("vec", vec)
	}
	scalar selc`i'=sqrt(vlc[1,1])
	scalar seec`i'=sqrt(vec[1,1])


	
	drop d_i
}

}
end


simulate estlc25=estlc25 estlc30=estlc30 estlc35=estlc35 estec25=estec25 estec30=estec30 estec35=estec35 selc25=selc25 selc30=selc30 selc35=selc35 seec25=seec25 seec30=seec30 seec35=seec35 true_25=true_25 true_30=true_30 true_35=true_35 , reps(1000): SPR

gen bias_estlc25=abs(estlc25-true_25)
gen bias_estlc30=abs(estlc30-true_30)
gen bias_estlc35=abs(estlc35-true_35)

gen bias_estec25=abs(estec25-true_25)
gen bias_estec30=abs(estec30-true_30)
gen bias_estec35=abs(estec35-true_35)

sum bias_estlc25 bias_estlc30 bias_estlc35
sum bias_estec25 bias_estec30 bias_estec35

sum estlc25 estlc30 estlc35 
sum estec25 estec30 estec35

sum selc25 selc30 selc35
sum seec25 seec30 seec35

foreach x in estlc estec{
	forv i=25(5)35 {
		qui sum bias_`x'`i'
		scalar bias=r(mean)
		qui sum `x'`i'
		scalar sd=r(sd)
		scalar rmse_`x'`i'=sqrt(bias^2+sd^2)
	}
}

disp rmse_estlc25
disp rmse_estlc30
disp rmse_estlc35

disp rmse_estec25
disp rmse_estec30
disp rmse_estec35

save "F:\econometrics\WP41_RatioNormEff\simulation\simulation_results_n40000.dta", replace











